Null infinity waveforms from extreme-mass-ratio inspirals in Kerr spacetime 
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We describe the hyperboloidal compactification for Teukolsky equations in Kerr spacetime. We 
include null infinity on the numerical grid by attaching a hyperboloidal layer to a compact domain 
surrounding the rotating black hole and the orbit of an inspiralling point particle. This technique 
allows us to study, for the first time, gravitational waveforms from large- and extreme-mass-ratio 
inspirals in Kerr spacetime extracted at null infinity. Tests and comparisons of our results with 
previous calculations establish the accuracy and efficiency of the hyperboloidal layer method. 
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I. INTRODUCTION 

Black hole perturbation theory plays a prominent role 
in obtaining physical insight into the quantitative behav- 
ior of astrophysical systems In recent years there 
has been considerable interest in computing the grav- 
itational waveforms emitted by the radiation-reaction 
driven inspiral of a small, compact object into a large 
supermassive black hole. The waveforms from such 
extremc-mass-ratio inspirals (EMRIs) are expected to be 
observed by future, space-based, gravitational wave de- 
tectors such as the Laser Interferometric Space Antenna. 

EMRIs provide a rich phenomenology relevant for as- 
trophysical and theoretical questions. Observations of 
EMRIs would deliver, among others, an accurate map- 
ping of black hole spacetimes @ and tests of the black 
hole uniqueness theorem [T^]. On the theoretical side, 
modeling of EMRIs poses challenging problems. Because 
the numerical solution of full Einstein equations for ex- 
treme mass ratios of the order of 10~^ are computation- 
ally prohibitive, EMRIs are calculated via approximation 
techniques. Currently, the most common approximative 
methods to study EMRIs are the kludge waveform ap- 
proach [nHIll, the effective-one-body (EOB) approach 
[H-Iil, and the self-force approach 

In these methods, the waveforms can be calculated by 
numerically solving the equations describing black hole 
perturbations. When the central, supermassive black 
hole is assumed to be rotating, the most common ap- 
proach is to solve the Teukolsky equations that describe 
curvature perturbations of a Kerr spacetime 0, [H] . 

A difficulty in solving the Teukolsky equations 
numerically — especially in the time domain — is the con- 
struction of suitable outer boundary conditions. Typi- 
cally, the computational domain is artificially truncated 
to calculate the solution in a finite domain. One needs 
to provide boundary data along the timelike boundary 
of this domain such that the artificial boundary is trans- 
parent to outgoing waves from the interior. The outer 
boundary problem is already nontrivial for the evolution 
of scalar waves in Minkowski spacetime j3ll - [33| . There 
has been substantial progress in constructing good outer 
boundary conditions for wave equations in Schwarzschild 
spacetime [33 - l37| . which have been recently generalized 
to Kerr spacetime [ssj . Such sophisticated boundary con- 



ditions, however, can be difficult to implement. To our 
knowledge, there is no numerical implementation of the 
improved boundary conditions in Kerr spacetime. 

Another difficulty in the numerical solution of the 
Teukolsky equation is the extraction of gravitational 
waves as measured by idealized observers at null infin- 
ity. Astrophysical sources of gravitational radiation are 
typically thousands of light years away, whereas the com- 
putational truncation is performed at a few thousand 
Schwarzschild radii. One argues that the dynamics in 
the asymptotic domain is negligible. However, there are 
certain effects, such as the backscattering off curvature 
jsol l40j , where asymptotic properties of the solution are 
essentially different from the corresponding properties at 
finite distances. Such differences can be relevant to gravi- 
tational wave astronomy. For example, it has been found 
that the asymptotic formula relating the curvature per- 
turbation -04 to the gravitational wave strain (Eq. (j25p ) is 
invalid for polynomially decaying solutions even at large 
distances used in numerical wave extraction [4l| . Fur- 
ther, we know that the phase of the gravitational wave 
signal measured at null infinity and at finite distances 
may difl:er substantially in EMRIs [H,!!!. 

A clean solution to both the outer boundary and the 
radiation extraction problems is to include null infinity in 
the computational domain. No artificial boundary condi- 
tions are needed when the solution is calculated globally 
in space, and the idealized gravitational wave signal can 
be read off directly at null infinity. 

There are two methods to include null infinity in the 
computational domain: the characteristic method and 
the hyperboloidal method. These methods correspond 
to approaching null infinity along null and spacelike di- 
rections respectively. 

The characteristic method is very well developed [3l ■ 
It has been applied recently to the unambiguous com- 
putation of gravitational waves from binary black hole 
mergers simulated via Einstein equations p5l - |47| . Con- 
cerning the Teukolsky equation, it has been used in 
Schwarzschild spacetime for black hole collisions in the 
close limit [3 IS| • The construction of a null foliation 
in Kerr spacetime, however, is rather complicated [sol - 
[53| . As a result, there is no numerical computation of 
gravitational perturbations of Kerr spacetime at null in- 
finity. 
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The second approach to including null infinity in the 
computational domain is the hyperboloidal approach 
[s^ ]. The construction of suitable hyperboloidal folia- 
tions is simpler than characteristic ones because spacelike 
surfaces are more flexible than null surfaces [s^ . Various 
studies have demonstrated the accuracy and efficiency 
of the hyperboloidal method in dealing with perturba- 
tions of black hole spacetimes, but most of them arc re- 
stricted to spherical symmetry. For example, numerical 
studies have been performed on Yang-Mills equations in 
a Schwarzschild spacetime (sgI \5T\. spherically symmetric 
wormholes [s^ and boson stars ISOf. self- force computa- 
tions [S^l and gravitational perturbations [4l|, |6l| includ- 
ing extreme-mass-ratio inspirals [H, li^l . 

The only numerical computations of black hole pertur- 
bations using the hyperboloidal method without spheri- 
cal symmetry deal with scalar perturbations [63 - l65| . One 
goal of this paper is to present the application of the 
hyperboloidal method to Teukolsky equations in Kerr 
spacetime for the calculation of gravitational waveforms 
at null infinity. The test problem we study to demon- 
strate the method is of astrophysical interest. The grav- 
itational waveforms we compute are emitted by the per- 
turbations due to a point particle inspiralling into the 
central Kerr black hole. 

Another goal of this paper is to solve a difficulty of ap- 
plying the hyperboloidal method in Kerr spacetime. The 
hyperboloidal foliation of Kerr spacetime constructed in 
[55| uses a transition zone in which a truncated Cauchy 
surface is matched to a hyperboloidal surface. This tech- 
nique requires tuning a large number of foliation parame- 
ters and results in a sudden drop of characteristic speeds 
in the transition zone decreasing computational efficiency 
[41I, IIO, HH . There are two solutions to this problem. One 
is to use a single smooth surface avoiding the transition 
zone of (ssf . This technique has recently been applied 
by Racz and Toth in a detailed study of polynomial de- 
cay rates of a scalar field in Kerr spacetime jg^]. They 
construct the first smooth, horizon-penetrating, hyper- 
boloidal foliation of Kerr spacetime for their study. Mod- 
ifying the coordinates everywhere smoothly leads to very 
efficient numerical computations. However, it may be fa- 
vorable due to practical reasons to employ standard coor- 
dinates in the strong field region near the rotating black 
hole, especially when particles arc present in the compu- 
tational domain, so that the description of the particles' 
motion docs not need to be changed. 

In this paper we use the recently proposed hyper- 
boloidal layer ([6^) instead of a transition zone. The new 
layer technique has been applied in studies of EMRIs us- 
ing the Regge-Wheeler-Zerilli formalism completed by 
EOB-resummed analytical radiation reaction (at linear 
order in the mass ratio) [i^, In this paper, we 

demonstrate the construction and application of a hy- 
perboloidal layer for the Teukolsky formalism and the 
kludge approach in Kerr spacetime. 

The paper is organized as follows. In Section II we 
present the standard Teukolsky formalism with the trans- 



formations typically used in numerical calculations. In 
Section III we present how the hyperboloidal layer tech- 
nique is applied to Kerr spacetime and to the Teukol- 
sky equations as two simple coordinate transformations. 
Numerical results obtained with this formalism are pre- 
sented in Section IV. After describing briefly the numer- 
ical method, we present a comparison of waveforms, an 
improved agreement of the energy fluxes with frequency- 
domain computations for several circular-equatorial or- 
bits, a discussion of recoil velocities, and an example sim- 
ulation that demonstrates remarkable gains in efficiency 
with the new method. A discussion of our results and an 
outlook can be found in Section V. 



II. THE TEUKOLSKY FORMALISM 

Curvature perturbations of Kerr spacetime are gov- 
erned by a separable equation in the Newman-Penrose 
formalism 0, Q . The original calculations involve the 
Boycr-Lindquist representation of the Kerr metric. It is 
convenient for numerical applications to introduce cer- 
tain transformations of the Teukolsky equation. We re- 
view the common transformations to set up the equations 
to which the hyperboloidal method will be applied. 



A. Teukolsky equation in Boyer— Lindquist 
coordinates 

The Kerr metric in Boyer-Lindquist coordinates 
{t, r, 9, if) reads 
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/ 2Mr\ , , AaMr n „ , , E , „ 
I dt^ sin^ edtdw+ — dr^ 



-hS de^ + r 



2Ma'^r sin^ 9 



sm^9dip^, (1) 



where E 



1^ cos 6*2, and A := + - 2Mr. We 



denote the mass of the Kerr spacetime by M, its angular 
momentum by aM. 

The Teukolsky equation describes curvature perturba- 
tions with spin weight s in an adapted Newman-Penrose 
orthonormal frame. The homogeneous Teukolsky equa- 
tion in Boyer-Lindquist coordinates reads 



A ^ 

■deismBde'i') + 



1 ^^-!^]d^^ 



_ f^M(c^_r^ + (r + m cos 0) ) 
2 ,a(r-M)^ zcos^. 



A sin^ 

- {s^ co\? 9 - s)^ , 

where D'^ = (j-^ + a^)^/A — a^sin^ 9. In the following, 
we review the transformations that are typically applied 
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to the Teukolsky equation for numerical computations 




B. Transformations of the Teukolsky equation 

Boyer-Lindquist coordinates have undesirable features 
for numerical computations. Their time slices intersect 
at the bifurcation sphere leading to a singular radial co- 
ordinate at the horizon, in a similar way as Schwarzschild 
time slices do. There are two ways to deal with the coor- 
dinate singularity at the horizon: transform to a horizon- 
penetrating time slicing, or push the horizon to infinite 
coordinate distance with the tortoise coordinate. 

Horizon-penetrating coordinates are regular at the 
event horizon. They are typically used in combination 
with excision. The numerical treatment of the inner 
boundary is clean with this method because there are 
no characteristics coming out of the horizon. In addi- 
tion, the radiation absorbed by the black hole can be 
calculated accurately. 

The most common approach in perturbation theory, 
however, is to use the tortoise coordinate. The tortoise 
coordinate pushes the horizon to infinite coordinate dis- 
tance. The infinite domain is then truncated at a finite 
distance, where artificial ingoing boundary conditions are 
imposed. These lead to a contamination of the physi- 
cal solution due to numerical reflections from the inner 
boundary. But these reflections are small because the 
potential terms in the Teukolsky equation fall off expo- 
nentially fast in the tortoise coordinate, so the ingoing 
boundary conditions are quite accurate. In our study, 
we focus on the outgoing radiation and accept small re- 
flections from the inner boundary. The truncation of the 
inflnite domain at a finite coordinate implies that the ef- 
ficiency of the code is not optimal because a large region 
in negative tortoise coordinate needs to be included in 
the computational domain. 

The relation of the tortoise coordinate to horizon- 
penetrating coordinates near the event horizon is sim- 
ilar to the relation of standard Cauchy coordinates to 



hyperboloidal coordinates near null infinity [4l|, A 
fundamental difference, however, is that the structure of 
the solution in the strong field region can be interesting, 
whereas the asymptotic solution is essentially simple in 
asymptotically flat spacetimes. The additional resolu- 
tion that the tortoise coordinate provides near the black 
hole can be relevant in certain studies. This question 
should be analyzed in more detail. We choose to employ 
the tortoise coordinate, and leave the study of horizon- 
penetrating coordinates to future work. 

The angular coordinate also needs to be transformed 
to cure the infinite twisting near the horizon |63| • Fur- 
ther, a rescaling respecting the fall-off behavior of the 
curvature perturbations is necessary. Finally, we sepa- 
rate the azimuthal dependence to arrive at a 2 -I- 1 di- 
mensional system. These transformations are listed as 
follows: 

1. Introduce the tortoise coordinate r*, 

dr ~ A ■ ^ ' 

2. Introduce the azimuthal coordinate 

dip = dip + -^dr. 

3. Rescale the curvature perturbation of spin weight 
s according to its asymptotic fall-off behavior, 

4. Separate the azimuthal dependence using the mode 
number k, 

In the following, we drop the subscript k from V-'fe 
for conciseness of notation. The Teukolsky equation be- 
comes [ill 



lA sm y 

2 

) + rsA + i{asA cos 9 + 2Mark)) dtip + 

+ 4r (A(8r2 + Qa^) + 2rs{r'^ + a'^){r - M) + 2iakr(r^ + a^)) dr ip + 
rA ■ ' 
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These are standard transformations for numerical com- 



putations. Hyperboloidal compactification adds only two 
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additional transformations to this list: a transformation 
of the time coordinate, and a compactification of the ra- 
dial coordinate. 



III. A HYPERBOLOIDAL LAYER FOR THE 
TEUKOLSKY EQUATION 

We present the hyperboloidal layer method in the lan- 
guage of standard transformations to make its applica- 
tion as straightforward as possible for future studies. 



A. Hyperboloidal compactification 

A hyperboloidal surface is an e very where spacelike sur- 
face that approaches null infinity 54| . It has been known 
already in the early days of numerical relativity that such 
surfaces are beneficial for studying gravitational radia- 
tion [toMts} . a systematic study of the hyperboloidal 
initial value problem started with Friedrich's work show- 
ing that future null infinity is smooth for certain nontriv- 
ial classes of initial data 7J] . There has been substantial 
effort to numerically simulate such data using Friedrich's 
conformally regular field equations [tsI [76j . 

Even though the benefits of hyperboloidal surfaces for 
radiation problems beyond the Einstein equations have 
been clear, it took more than two decades before the 
hyperboloidal method could be applied successfully to 
systems other than the conformally regular field equa- 
tions. The first application of this kind is the study of 
magnetic monopoles by Fodor and Racz [tJ [iS ■ Their 
hyperboloidal foliation is based on a scri-fixing gauge in 
Minkowski spacetime first explicitly constructed by Mon- 
crief [ll] (see also [Z^). It was suggested that scri-fixing 
gauges should be beneficial also for studies on black hole 
spacetimes, however, the construction of a good hyper- 
boloidal coordinates turned out to be difficult |8C)I - [85| . 

The general construction of suitable hyperboloidal, 
scri-fixing gauges has been presented in [5q . We follow 
[ssf to present the hyperboloidal compactification. The 
technique consists of introducing a compactifying radial 
coordinate p and a suitable time coordinate r. 



The radial compactification can be performed conve- 
niently with the following definition of a compactifying 
coordinate: 

(4) 



1. Radial compactification 

We map the infinite positive domain in the tortoise 
coordinate r* to a finite domain in the compactifying 
coordinate p to compute the solution in an unbounded 
physical domain using a finite numerical grid. After com- 
pactification, the outer boundary of the numerical grid 
corresponds to infinity with respect to the physical coor- 
dinate. Infinity has a rich structure in spacetime. What 
part of infinity is included in the computational domain 
with radial compactification depends on the nature of the 
time surfaces, as discussed in the next subsection. 



The choice of Vl determines the properties of the com- 
pactification. Its zero set corresponds to infinity with re- 
spect to r*. For example, a common choice encountered 
in the literature is = 1 — p. We need more freedom in 
the choice of compactification, therefore we will state the 
general properties that the function f2 needs to satisfy. 
These properties are those that are satisfied by a con- 
formal factor in the Penrose conformal compactification 
picture [IB, 1131 ■ For all finite we require O > 0. De- 
noting the zero set of as S" (the p-coordinate location 
of r*-infinity), we require 

0(5') = 0, n'{S)^Q. 

where fi' = dfl/dp. 

It is known that compactifying the radial coordinate 
along Cauchy surfaces results in resolution problems [s^ . 
These problems can be avoided by a hyperboloidal time 
transformation that changes the compactification at spa- 
tial infinity to a compactification at null infinity. 



2. Time transformation 

A suitable time transformation shall leave the exterior 
timelike Killing field invariant so that no gauge dynam- 
ics is introduced into the stationary background. The 
new time coordinate r satisfies then dr = dt, which is 
achieved by a transformation of the form t = t — h{x'^), 
where are the spatial coordinates. The function h is 
called the height function. Most useful foliations of sta- 
tionary spacetimes are given via such a relation. We are 
interested in radially outgoing gravitational waves, there- 
fore we choose the height function to depend only on the 
tortoise coordinate 



T = t — h{r^,). 



(5) 



The height function is chosen such that the surfaces r = 
const, are hyperboloidal [s^. The specific choice we use 
is presented in the next section (see Eq. ([7])). 

In summary, hyperboloidal compactification consists 
of two simple transformations that can be written in the 
form ^ and ^ with free functions fl and h. The specific 
choice of these free functions depends on the background 
spacetime and its coordinatization. There is a large free- 
dom in their choice which allows us to use standard co- 
ordinates in an interior domain as described below. 



B. Hyperboloidal layer for Kerr spacetime 

An essential advantage of the hyperboloidal method 
is that it requires modifications only in the asymptotic 
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domain, which allows us to use arbitrary coordinates in 
the strong field regime. 

A transition zone is used in [ssf for matching arbitrary 
coordinates near the central black hole to hyperboloidal 
scri-fixing coordinates in the asymptotic domain. Subse- 
quent numerical experiments showed that the transition 
zone requires the adjustment of many arbitrary parame- 
ters and leads to a sudden drop of outgoing characteris- 
tic speeds decreasing the efficiency of the hyperboloidal 
method [ill, [e^, [sl . 

A new technique that avoids these deficiencies is the 
hyperboloidal layer presented in [6^. The idea is to at- 
tach a radial shell to the truncated computational domain 
in which the foliation is hyperboloidal and the radial co- 
ordinate is compactified. The layer needs to be attached 
in a sufficiently smooth way for the stability of the nu- 
merical computation. 

We would like to have the new coordinate p agree with 
the tortoise coordinate in the strong field region. De- 
noting the location of the interface between the interior 
domain and the layer by i?*, we set J7 = 1 for p < R^, 
so that p = in the interior. We also require that the 
coordinates agree up to a certain order at the interface. 
These conditions arc fulfilled by the following expression 



f7 = 1 



p__R* 



e(p-i?,): 



(6) 



where 8 denotes the step function and S denotes the co- 
ordinate location of the outer boundary. We emphasize 
that various similar choices for the radial compactifica- 
tion are possible. 

The original construction of a height function for a 
hyperboloidal layer [66j requires unit characteristic speed 
in the layer with respect to the layer coordinates r and 
p. An equivalent but more intuitive requirement is that 
outgoing null rays have the same representation in the 
layer coordinates as in the interior coordinates [43| . In 
Schwarzschild spacetime with the tortoise coordinate the 
requirement reads t — p = t — r*. From this relation 
together with ^ we get 



nip) 



(7) 



In the domain where = 1, we have /i = so that we 
obtain the standard coordinates. For the derivative of 
the height function we get 

dr^, dr^, ft — pil' ^ ^ 

We refer to H as the boost function. It vanishes in the 
interior domain where standard coordinates are used with 
r2 = 1; it is unity at infinity, where = 0. 

The Refs. [H, use the prescription ^ in 

Minkowski spacetime with standard coordinates, and in 
Schwarzschild spacetime with the tortoise coordinate. In 



this paper we apply it for the Kerr spacetime in Boyer- 
Lindquist tortoise coordinates. This procedure leads to 
a regular hyperboloidal compactification in Kerr space- 
time, even though the outgoing null rays do not have 
the simple form t — r^,, because the Kerr metric in the 
tortoise Boyer-Lindquist coordinates has asymptotically 
the same form as the Schwarzschild metric in tortoise 
Schwarzschild coordinates, or the Minkowski metric in 
standard coordinates. This feature is another evidence 
for the fiexibility of hyperboloidal coordinates. Only the 
asymptotic behavior of the null cone is relevant for the 
regularity of hyperboloidal compactification, whereas in 
the characteristic method the exact local form of null 
cone plays an essential role, thereby restricting the coor- 
dinate transformations [SOl - ISSj . 

The choices Q and ^ with the free functions explic- 
itly given in (jG]) and ([7]) completely fix the coordinate 
transformation from {t,?"*} to {t, p} that gives us a hy- 
perboloidal layer in Kerr spacetime for the coordinate 
domain [R^ , S] . 



C. Hyperboloidal compactification of the 
Teukolsky equation 

The transformations discussed in the previous section 
can be regarded as additional items in the list of trans- 
formations for the numerical solution of the Teukolsky 
equation given in Section III Bl 

5. Introduce a compactifying coordinate p, 

P 



6. Introduce a time coordinate t, 

T = t — h{r^). 



(9) 



(10) 



The compactification function fl{p) and the height func- 
tion /i(r*) are determined via (|6]) and ([7]). With the co- 
ordinate transformations at hand, we can now proceed 
to transform the Teukolsky equation. 

It is a priori not clear that a simple coordinate trans- 
formation will lead to a regular hyperboloidal compact- 
ification of the Teukolsky equations. There are various 
possibilities to study the regularity of such compactifica- 
tion. One is to derive conformal Teukolsky equations 
with respect to a conformally extended, regular Kerr 
metric. It has been shown in a previous application of 
hyperboloidal compactification using the Teukolsky for- 
malism in Schwarzschild spacetime (called the Bardeen- 
Press equation after that such an analysis is not nec- 
essary [61|. However, the method of [6l[ still requires 
the calculation of the Teukolsky equations ab initio in 
a general orthonormal Newman-Penrose frame which is 
then adapted to a hyperboloidal scri-fixing gauge to en- 
sure regularity of hyperboloidal compactification. The 
method we use in this paper is easier to apply. 
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A general operator for an equation in the independent 
variables {t,r^,6}, such as the Teukolsky equation 
can be written as 0[ip] — 0, where O is the differential 
operator 

O = A"d^ + A'^-'dtdr, + A'-'^-^dl^ + A^^dl 

The cocfhcients depend on r^, and 9] they can be read 
off from ([3]). With the transformations © and ((TU)) the 
derivative operators with respect to the old coordinates 
can be written in terms of the new coordinates as 



dt = dr. 



dr 



-Hdr + {1- H)dp. 



We used ([5]) in the expression for d,., ■ We write the 
transformed operator as 



O = A^^d' + A'^Pdrdp + APPdi + A'^^df 



-B^dr + BPdp + 



(11) 



where 



j^rr 


= A" 


- HA'''' - 


hA'-'^'-H^, 




= (1- 


-H){A''-' 


-2H A' 


APP 


= (1- 


- H fA''''' 




B^ 


= B* - 


- HB''' - 


A'-'''H'{1-H), 


BP 


= (1- 


-H)iB'-'- 


-A^'-'^'H'), 



and H' := dH/dp. All coefficients are functions of p, and 

e. 

The transformed system is by construction equiva- 
lent to the original system ([3]) in the interior domain 
p < i?* where we have H ^ 0. We need to make 
sure that the hyperboloidal compactification is regular 
at infinity where H = 1. Note that asymptotically 



(l-H) 



We see by inspection that 



the coefficients of the Teukolsky equation ([3]) have the 
correct asymptotic behavior ensuring regularity of the 
hyperboloidal compactification. The transformed coeffi- 
ciencts have explicitly finite values at future null infinity. 
We mention that this feature is not special to the Teukol- 
sky equations. Similar regular hyperboloidal compacti- 
fications for other wave equations have been studied as 
mentioned in the Introduction. 

We also need to make sure that the transformed system 
is pure outflow at the outer boundary so that no bound- 
ary conditions are required. This condition, along with 
the regularity, can be checked by evaluating the trans- 
formed equations at infinity, that is, at {p = 5}. We 
obtain for s = — 2 



-S(S - R,)drdpij + 2dli) 

~4:{UI + ia{k - 2 cos e))dr^p 

+2 cot edetp - 2(3 + fc^ - 4fc cos 6 + cos(26'))V', 



where D = S{S — i?*) — 2a^ sin^ 9. The principal part of 
the operator is of the form Ci dridr + dp) + C2 d"^ + 2dg, 



where Ci and C2 are coefficients. Therefore, the modes 
of the system propagate either along the boundary at 
{p = S}, or out of the domain along the level sets of 
characteristics t — p. As expected, there are no incoming 
modes. 



IV. NUMERICAL RESULTS 

In this section, we present tests using the hyperboloidal 
compactification of the Teukolsky equation. Our main 
result is the efficiency and accuracy of the hyperboloidal 
layer method providing unambiguous waveforms at fu- 
ture null infinity from large- and extreme-mass-ratio in- 
spirals in Kerr spacctime. 



A. The numerical method 

For the discretization of the continuous equation we 
use a two-step Lax-Wendroff algorithm as in [63] ■ After 
performing the transformations presented in the previous 
sections, we rewrite the vacuum equation in the form 

d'^tP = A'-PdrdptP + APPd^ptP + A'^'^d^^p 

+ B^dryJ + BPdp^p + B^dgyj + C^p, (12) 

where the coefficients with a tilde are obtained by divid- 
ing the coefficients of the operator in pT|) by —A'^'^. We 
put the equation (jl2p in first-order form (in p and r) by 
defining a new field variable 



TT = drip + bdpip , 



b = -{A^P + J{A^P)2 +4APP)/2 



(13) 
(14) 



We chose to use this first-order form of the equation for 
our numerical evolutions, because it has been discovered 
through extensive experimentation [g^I that such a form 
is ideally suited for long stable evolutions. Then equation 
(IT2]) with source terms T takes the form 



drU 



MdpU 



Lu + Au 



where 



U : 



{tP'R,1pI,Trji,TTl} 



(15) 



(16) 



is the solution vector. The subscripts R and I refer to 
the real and imaginary parts respectively. The matrices 
M, A and L are obtained from (fT2|) as in [l^. Here 
it will suffice to simply indicate the final form taken by 
these matrices: 



M = 



f b \ 

6 

"131 -b 

-77132 77731 -bj 



(17) 
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0-1 


131 0^32 ^33 



(18) 



V- 



and 



032 031 —034 033^ 



o^ 
0000 
^31 000 

^31 Oy 



(19) 



The angular derivatives are encoded in L. 

The main advantage of casting the equation (|12p into 
the form psp is that the system has advantageous proper- 
ties in the variable p. The matrix M has a complete set of 
hncarly independent eigenvectors with real eigenvalues. 
This is not a rigorous statement on the hyperbolicity of 
the system because the matrix L contains second-order 
angular derivatives. Nevertheless, experiments show that 
the system is numerically well-behaved. A time-explicit 
evolution scheme is developed for this system of linear 
partial differential equations using the two-step, second- 
order Lax-Wcndroff finite-difference method. We write 
(fist as 



drU + DdpU = S , 



where 



D 



/6 

6 

0-60 

Vo -b/ 



S = -(M - D)dpU Lu Au + T. 



(20) 



(21) 



(22) 



Each iteration consists of two steps. In the first step, the 
solution vector between grid points is obtained from 
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This is used to compute the solution vector at the next 
time step. 
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The angular subscripts arc dropped in the above equation 
for conciseness. All angular derivatives are computed us- 
ing second-order, centered finite difference expressions. 
Symmetries of the spheroidal harmonics are used to de- 
termine the angular boundary conditions: For even \m\ 
modes we have deijj = at = 0, tt; for odd \m\ modes we 
have V' = at 6* = 0, TT. We impose "ingoing" boundary 
conditions at the inner radial boundary. We do not need 
to apply any boundary condition at the outer boundary. 



The small compact object inspiraling into the central 
supermassive black hole is modeled as a point particle in 
the large- or extreme- mass-ratio limit. The representa- 
tion of point particle in mathematics is performed typi- 
cally by a Dirac delta distribution. More specifically, the 
particle energy-momentum tensor takes the form 



Sf sin 6 



S[r, ~ r,{T)]S[e - e{T)]S[^ - cj>{T)] 



where Uq and p refer to the 4-vclocity and the rest mass 
of the particle, respectively. Note that f = dr/ds ap- 
pears in the denominator of this expression, which tends 
to 00 as the particle approaches the horizon (we denote 
by s the proper time along the particle's trajectory). 
Thus, the particle source-term on the right-hand-side of 
smoothly decays away as the particle approaches 
the horizon, thereby allowing the evolution equation to 
gradually transition into its homogeneous form. This 
source-term is constructed using the energy-momentum 
tensor describing a point particle moving in Kerr space- 
time depicted above. The explicit expression for T is 
very lengthy and not particularly illuminating. Here, we 
simply point out that the final expression is built us- 
ing Dirac delta functions in r and 9, as well as first and 
second derivatives of the delta functions in these vari- 
ables. These terms have coefficients that are complex 
functions of the black hole's physical parameters and also 
the trajectory of the point particle. Details on the par- 
ticle source-term, along with the representation of the 
delta function and its derivatives on a numerical grid, 
are given in [iMl,[8|. 

The trajectory of the particle in a decaying orbit 
around the central black hole is computed separately 
and is then used to calculate the source-term mentioned 
above, which in turn is fed into the Teukolsky equation 
solver code. The computation of the decaying trajectory 
can be broken into three distinct pieces: An early time 
adiabatic inspiral, in which the particle is approximated 
as evolving through a sequence of bound orbits of the 
central black hole that are calculated using an "energy 
balance" approach; a late time geodesic plunge, in which 
the particle falls into the black hole and radiation re- 
action is ignored; and an intermediate transition regime 
that smoothly connects these two pieces. Details on how 
these steps are handled and use d by the Teukolsky equa- 
tion solver code can be found in [l5| . It is straightforward 
■ to make use of decaying orbital trajectories from other 
approaches, such as the EOB or the self-consistent self- 
force approaches. 



B. Waveforms 

After solving the Teukolsky equation with a parti- 
cle source in the time-domain, we compute the gravita- 
tional waveforms and using the simple relationship 



8 



shown below that is vahd in the far field, 



9^ 



(25) 



Note that in our study we have direct access to the far 
field because null infinity is part of the computational 
domain. This allows us to use the above relation as an 
equality and to extract the waveform cleanly. 

Figures [T] and [5] depict the gravitational waveform him 
measured at future null infinity, with £ ~ m ~ 2. The 
mass-ratio for this evolved binary system is /i = 10^'^ and 
the spin of the central Kerr black hole is a = 0.5. The sys- 
tem undergoes a decay due to gravitational wave emission 
over approximately 123 full cycles. The computational 
domain is p e [—50,50] and 9 E [0,7r] with grid sizes of 
3125 and 32 respectively. The hyperboloidal layer starts 
at yO = 14, therefore the particle's orbit never crosses into 
the layer. There arc a total number of 562, 500 time-steps 
in this computation. 

The figures depict the waveforms from both the orig- 
inal Cauchy code [ij, and the new code with the 
hyperboloidal layer in Kerr spacetime developed in this 
work. The Cauchy code approaches spatial infinity, and 
therefore cannot provide direct access to waveforms at 
null infinity. To obtain data from the Cauchy code to 
compare with that from the hyperboloidal layer code, we 
extract the waveforms at multiple radii on the spatial 
grid and extrapolate them along fixed values of retarded 
time, to infinity, using a simple second order polynomial 
expansion. We do not employ higher order extrapola- 
tion for the Cauchy code because the truncation error 
becomes larger than extrapolation error with the given 
resolution and second order finite differencing. Never- 
theless, the waveforms from the two codes agree so well 
that it is difficult to make out any difference between the 
two in these figures. In Figure [3] we show the differences 
in the waveform amplitude and phase for the same data. 
The relative difference between the amplitudes is at the 
level of 10"'^ and the maximum difference in phase is 
0.06 rad. A better agreement between extrapolated and 
null infinity waveforms can be obtained with more accu- 
rate extrapolation algorithms and higher-order methods 
in spherical symmetry j43| . 



C. Fluxes 

Gravitational waves carry energy away from a binary 
system thus causing it to inspiral. The calculation of the 
energy carried away by these waves uses the standard 
expression, 



dE 1 
— — = lim — 

dr r-i-oo 4:77 



^dt' 



dfl 



(26) 



in a postprocessing step after time evolution. In standard 
computations, the limit to infinite radius in the above 
equation is approximated by application of the formula 
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Figure 1: The full waveform amplitude emitted by the in- 
spiral, plunge, and ringdown of the particle into the central 
Kerr black hole calculated both with the hyperboloidal code 
(black), and the Cauchy code (dashed, red). The merger time 
(r — 0) is defined rather arbitrarily by the time of maximum 
amplitude. 
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Figure 2: The waveform emitted during the final plunge phase 
of the EMRI inspiral with the largest gravitational radiation 
emission. 



at finite but large radii, and a subsequent extrapolation 
of the corresponding values to infinity. One major ad- 
vantage of the hyperboloidal method is that the limit in 
Eq. (pS)) can be replaced by local evaluation of the in- 
tegral along future null infinity, which is the boundary 
of the computational domain. Therefore, the extrapo- 
lation of fiuxes is not necessary, and the energy can be 
computed cleanly. There are more subtleties in the com- 
putation of the fiuxes related to the integration in time 
from the infinite past We ignore these and set the 
integration constant to zero. We observe empirically that 
this procedure does not lead to large errors. 

In Table |I] we show numerical values of the energy 
fluxes computed at null infinity by the new hyperboloidal 
layer code for several circular-equatorial orbits. Com- 
paring these values with those from very high-accuracy 
frequency-domain approaches [9l| , we note an agreement 
at the 99.95% level. We interpret this high level of agree- 
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into that layer throughout the simulations. 



Figure 3: Relative difference in amplitude (top) and phase 
difference (bottom) between the hyperboloidal code and the 
Cauchy code. For clarity, it is worth pointing out that the 
phase difference reaches a maximum in the ballpark of 0.06 
rad near r = 0. 



a 


Radius 


TD 


FD 


% Diff. 


Conv. 


-0.5 


10.576 


4.5560 


4.5584 


0.053 


1.96 


+0.0 


6.0000 


73.683 


73.720 


0.050 


2.05 


+0.5 


4.2753 


284.72 


284.88 


0.056 


2.04 


+0.8 


11.627 


2.2272 


2.2287 


0.067 


2.04 



TABLE I: This table depicts the numerical values of energy 
fluxes for the \m\ — 2 mode from various circular-equatorial 
orbits, computed by the time-domain hyperboloidal code 
(TD) and a high-accuracy frequency-domain (FD) code. The 
flux values are scaled up with the square of the inverse mass- 
ratio (/i~^) and by a factor of 10^. Convergence rate is also 
presented. 



ment with the frequency-domain fluxes as strong evi- 
dence for the accuracy of the code, especially as com- 
pared with the original time-domain Cauchy code that 
achieves agreement with frequency-domain fluxes at the 
level of 99% .13, 8^. 

In Table U we also show the standard three level con- 
vergence rates. We obtain clear second-order conver- 
gence, as expected. The fluxes depicted in the table were 
obtained via Richardson extrapolation using the data 
from the different grid resolutions. The computational 
domain is p G [—50,50] and 6 E [0,7r] with grid sizes 
1250 X 32, 2500 x 64 and 5000 x 128. The hyperboloidal 
layer starts at p = 14, therefore these orbits never cross 




Figure 4: Linear momentum flux measured at {p = S} with 
the hyperboloidal code (black), and extrapolated to infinity 
with the Cauchy code (dashed, red). 
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Figure 5: Linear momentum flux as in Fig. |4]zoomed into the 
flnal plunge stage of the inspiral. 



D. Kick velocities 

The asymmetric emission of gravitational waves during 
the inspiral of the small black hole into the supermassive 
central black hole leads to a recoil due to conservation of 
momentum. We compute this recoil by first computing 
the linear momentum flux carried away by the gravitation 
waves^and then performing a simple time-integral of this 
flux [l5|. To compute the momentum flux, we use the 
expression 



lim - — 

7'— ^OO 4:71 



%l)dt' 



(27) 



where the is the unit direction vector in spherical co- 
ordinates. As for Eq. ((26|). the limit to infinity is replaced 
by local evaluation of the integral at future null infinity, 
which allows us to read off the momentum flux without 
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Figure 6; Linear momentum in x direction. 

an approximation in the radial direction. The integration 
constant resulting from starting the integration at r = 
is set arbitrarily to zero. 

In Figures |31 [S] and H] we depict the results from 
the computation of linear momentum flux and its time- 
integral. The mass-ratio for this computation is fi — 10~^ 
with the central black hole spin a = 0.7. The computa- 
tional domain is p € [—50,50] and 9 e [0,7r] with grid 
size 3125 x 32. The hyperboloidal layer starts at p = 14, 
therefore the particle's orbit never crosses into that layer 
throughout the simulation. There are 2, 250, 000 time- 
steps in this computation. Once again we superimpose 
results from the Cauchy code and the new hyperboloidal 
code which agree to a very high level. The final values of 
the recoil speed from the two different codes, 1.645 x 10~^ 
and 1.641 x 10~^, agree to three significant digits, con- 
sistent with our previous results. 

E. Efficiency 

Previous sections provided evidence for the accuracy 
of the hyperboloidal method. The gain in accuracy in 
observables, such as the waveforms or energy fiuxes, is 
mainly due to the direct numerical access to asymptotic 
quantities at future null infinity, which is part of the com- 
putational domain. The difference in those quantities, 
however, is rather small, essentially because extrapola- 
tion in given background spacetimes already gives reli- 
able results. The basic advantage of the hyperboloidal 
method in this context is that it simplifies the accurate 
extraction of physical observables from numerical simu- 
lations: It eliminates the extrapolations. 

The remarkable feature of the hyperboloidal method 
is that it provides improved accuracy at negative cost. 
The numerical computation is much more efficient along 
hyperboloidal surfaces than along Cauchy surfaces. One 
may think that the calculation of waveforms at null infin- 
ity introduces additional computational expense for nu- 
merical simulations. In fact, the opposite is true: One 
gains more accuracy for less computational cost. 



The explanation of this property of hyperboloidal evo- 
lutions is that the standard method of using truncated 
Cauchy foliations with artificial outer boundary condi- 
tions is unnatural to study radiative solutions. An effi- 
cient computation of the asymptotic gravitational wave 
signal can only be achieved if the surfaces, along which 
the solution is computed, follow the outgoing signal 
closely. Gravitational waves propagate to infinity along 
null rays. Null rays and Cauchy surfaces diverge infinitely 
from each other in the asymptotic region. Cauchy sur- 
faces are therefore inadequate for computing outgoing 
asymptotic radiation. Further, their truncation at a fi- 
nite distance introduces errors through artificial outer 
boundary data contaminating the computation. 

To avoid such contamination in a standard Cauchy 
code the outer boundary is usually chosen to be causally 
disconnected from the region of interest. This implies 
that the duration of the simulation is limited by the size 
of the numerical grid. When null infinity is part of the 
computational domain, however, no such restriction ap- 
plies. Therefore, in principle, we can perform arbitrary 
long evolutions, only restricted by the accumulation of 
the numerical truncation and the physical approximation 
errors. The efficiency of a hyperboloidal code as com- 
pared to a Cauchy-type code increases with the duration 
of the simulation. 

To demonstrate the efficiency of our hyperboloidal 
Teukolsky code, we performed an evolution of an EMRI 
that lasts for one million M. The spin of the central 
black hole is 0.9, the mass ratio is 10~^. There are about 
10, 111 orbital cycles in this simulation. The computa- 
tional domain is p G [—50, 50] and 6 G [0, tt] with the grid 
size 5000 x 32. The hyperboloidal layer starts a.t p = 14, 
therefore the particle's orbit never crosses into that layer 
throughout the simulation. There are over a 100 million 
time-steps involved in this computation. In Fig. [7] we 
show the amplitude and frequency against time in a half- 
logarithmic scale. The amplitude and frequency of the 
gravitational waveform at null infinity is fairly constant 
over most of the simulation, as expected. The main con- 
tent of this plot is to demonstrate that such simulations 
in time-domain are now possible with the new method. 
Using the Cauchy code for such a simulation, would have 
required us to place the outer boundary at least as far as 
500, OOOM. So the Cauchy code would take 5000 times 
longer to perform this simulation with equivalent local 
grid resolution. This number can be increased even fur- 
ther using horizon-penetrating coordinates. 

A million M in geometric units for a supermassive 
black hole of, say, 6 million solar masses corresponds to 
about one Earth year in real time. The simulation for 
Fig. [7] has been performed on a computer cluster using 
on the order of 1000 processor-cores in just under a day. 
This shows that we can perform long-time simulations 
of EMRIs much faster than the expected real time du- 
ration of the events, and that large parameter studies 
in perturbation theory are computationally feasible with 
time-domain codes. 
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Figure 7: An evolution that runs a million M. The frequency 
is depicted in red. The amplitude (top black line) has been 
scaled up by The inset shows the final merger. 



V. DISCUSSION AND OUTLOOK 

We presented the first numerical evolution scheme for 
computing gravitational perturbations of a rotating black 
hole spacetime including future null infinity. We tested 
our method on the inspiral of a point particle into a cen- 
tral rotating black hole. The main advantages of the 
method are: (i) increase in accuracy by direct radiation 
extraction at future null infinity; (ii) a clean solution to 
the outer boundary problem in Kerr spacetime; and (iii) 
a great gain in efficiency for long time simulations. 

Our method is based on the hyperboloidal compact- 
ification of Kerr spacetime presented in [s^. We also 
showed, for the first time, how the hyperboloidal layer de- 
veloped in [66j can be employed in a Kerr spacetime. The 
layer technique allowes us to use arbitrary coordinates in 
an interior compact domain including the central black 
hole and the orbit of the particle, while compactifying the 
asymptotic domain along a hyperboloidal foliation. Be- 
cause the hyperboloidal layer is attached to the interior 
domain sufficiently smoothly, only minor modifications of 
existing computational infrastructure is needed (see [43| 
for a related study in Schwarzschild spacetime) 

We showed that the hyperboloidal compactification of 
the Teukolsky equations can be performed by a simple 
coordinate transformation. There is no need for study- 
ing a conformally regular Kerr spacetime, or for calcu- 
lating the equations ab initio in a general orthonormal 
Newman-Penrose frame as has been done in [6l|. The 
coordinate transformation of the equations as given in 
Boycr-Lindquist coordinates leads naturally to a regular 
hyperboloidal compactification. Only the asymptotic be- 
havior of the null cone is relevant for the hyperboloidal 
method, which implies that the simple prescription de- 
vised in Minkowski and Schwarzschild spacetimes can be 
used in Kerr spacetime with respect to Boyer-Lindquist 
tortoise coordinates. This provides evidence for the fiex- 
ibility of the hyperboloidal approach, as opposed to the 



characteristic approach which restricts the coordinates 
locally for compactification at future null infinity [50l - f53| . 

We compared our results to a previous study of EM- 
RIs in Kerr spacetime [l^. The emphasis in this pa- 
per is not on new physics but on the numerical accuracy 
and efficiency of the hyperboloidal method as applied in 
Kerr spacetime. The waveforms at null infinity differ only 
slightly from their previous computation obtained by ex- 
trapolation. We showed that the hyperboloidal method 
gives better accuracy in the energy fluxes at inflnity. Our 
strongest case for efficiency is the simulation of an inspiral 
that lasts one million M in geometric units, which would 
have taken about 5000 times longer using the standard 
method. 

In the future this method should be applied further 
to physically interesting problems. One could, for ex- 
ample, look for the new ringdown frequency mode for 
perturbations of a Kerr black hole discovered by Mino 
and Brink [ol] (see also [9^) or for transient resonances 
expected for certain orbits in Kerr spacetime |94| . Such 
studies require high accuracy, so the description of the 
particle's motion needs to be improved, for example by 
including conservative effects in the kludge approach, by 
reading off the source terms from an EOB description, or 
by computing the self-consistent evolution driven by the 
high-order in mass ratio self-force of the particle. 

It is clear that there is much room for development 
in the accurate computations of EMRIs. It would be 
interesting to compare the different approaches to the 
EMRI problem in the same setting, similar to compara- 
tive studies performed in nonlinear numerical relativity 
[95| . Comparison of physical invariants on given back- 
ground spacetimes is simpler, therefore such a study 
should be easier to perform in the perturbative setting 
than in the nonlinear one. 

On a numerical level, aside from the technical im- 
provements such as the use of high-order finite differenc- 
ing or multidomain pseudospectral methods, implemen- 
tation of horizon-penetrating coordinates near the black 
hole should lead to a further improvement in efficiency. 
Horizon-penetrating, hyperbolodial surfaces have certain 
advantages in perturbation theory over the standard 
Boyer-Lindquist or Schwarzschild time surfaces that in- 
tersect at the bifurcation sphere and approach spatial 
infinity (69j . In our context, horizon-penetrating coor- 
dinates should allow us to calculate accurately the ab- 
sorbed fluxes at the horizon via the ingoing radiation 
represented by the curvature component ip^. 

Our calculation of an inspiral that corresponds to 
roughly a year in real time (depending on the mass of the 
central supermassive black hole) but takes only about one 
day of computation is a strong evidence for the efficiency 
of our numerical infrastructure (see Section flV Ep . Note 
that the situation is reversed when fully nonlinear Ein- 
stein equations are solved for equal-mass, stellar, binary 
black-hole mergers. There, a typical computation takes 
about a month, but a stellar black-hole merger takes only 
milliseconds in real time. Therefore, it is imperative that 
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the efficiency of nonlinear Einstein codes is increased. 
Given the physical simplicity of an equal-mass binary 
black-hole configuration, a combination of technical im- 
provements including the hyperboloidal method, along 
with analytical insight into the problem, should allow us 
to simulate such systems much more efficiently in the 
near future. 

The application of the hyperboloidal method to generic 
computations with the fully nonlinear Einstein equations 
is an outstanding problem. Recently there has been stud- 
ies on various aspects of the hyperboloida l ap proach such 
as initial data |96|, or gauge conditions |97|. There are 
also suggestions on the evolution problem that do not 
require exp l icit regularity of the equations at future null 
infinity j98l - llOCl{ . The only successful numerical imple- 
mentation of such a formalism is by Rinnc in axisymme- 
try |lOl| , but no generic computation could be presented 
so far. We hope that the insight provided by studies in 
perturbation theory will help the application of the hy- 



perboloidal method to nonlinear Einstein equations. 
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